{ "cells": [ { "cell_type": "markdown", "id": "0", "metadata": {}, "source": [ "# Example: Extracting a state-space model from an equilibrium\n", "\n", "This notebook builds a self-contained MAST-U-like equilibrium, then constructs and inspects the FreeGSNKE state-space model. This is useful for magnetic control design." ] }, { "cell_type": "markdown", "id": "1", "metadata": {}, "source": [ "## Build the equilibrium" ] }, { "cell_type": "code", "execution_count": null, "id": "2", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import numpy as np\n", "\n", "# build machine\n", "from freegsnke import build_machine\n", "tokamak = build_machine.tokamak(\n", " active_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_active_coils.pickle\",\n", " passive_coils_path=f\"../machine_configs/MAST-U/MAST-U_like_passive_coils.pickle\",\n", " limiter_path=f\"../machine_configs/MAST-U/MAST-U_like_limiter.pickle\",\n", " wall_path=f\"../machine_configs/MAST-U/MAST-U_like_wall.pickle\",\n", ")\n", "\n", "# initialise equilibrium object\n", "from freegsnke import equilibrium_update\n", "eq = equilibrium_update.Equilibrium(\n", " tokamak=tokamak,\n", " Rmin=0.1, Rmax=2.0, # radial range\n", " Zmin=-2.2, Zmax=2.2, # vertical range\n", " nx=65, # number of grid points in the radial direction (needs to be of the form (2**n + 1) with n being an integer)\n", " ny=129, # number of grid points in the vertical direction (needs to be of the form (2**n + 1) with n being an integer)\n", " # psi=plasma_psi\n", ") \n", "\n", "# initialise profile object\n", "from freegsnke.jtor_update import ConstrainPaxisIp\n", "profiles = ConstrainPaxisIp(\n", " eq=eq,\n", " paxis=8.1e3,\n", " Ip=6.2e5,\n", " fvac=0.5,\n", " alpha_m=1.8,\n", " alpha_n=1.2\n", ")\n", "\n", "# initialise solver\n", "from freegsnke import GSstaticsolver\n", "GSStaticSolver = GSstaticsolver.NKGSsolver(eq) \n", "\n", "# set coil currents\n", "import pickle\n", "with open('data/simple_diverted_currents_PaxisIp.pk', 'rb') as f:\n", " current_values = pickle.load(f)\n", "\n", "for key in current_values.keys():\n", " eq.tokamak[key].current = current_values[key]\n", "eq.tokamak[\"P6\"].current += 500\n", "\n", "# carry out forward solve\n", "GSStaticSolver.solve(eq=eq, \n", " profiles=profiles, \n", " constrain=None, \n", " target_relative_tolerance=1e-9)\n", "\n", "# plot the resulting equilibrium\n", "fig1, ax1 = plt.subplots(1, 1, figsize=(4, 8), dpi=80)\n", "ax1.grid(True, which='both')\n", "eq.plot(axis=ax1, show=False)\n", "eq.tokamak.plot(axis=ax1, show=False)\n", "ax1.set_xlim(0.1, 2.15)\n", "ax1.set_ylim(-2.25, 2.25)\n", "plt.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Calculate state space model quantites\n", "\n", "Instantiate the nonlinear solver object in order to calculate the state space model matrices and growth rate parameters (and stability margin) for this equilibrium.\n", "\n", "Next you need to choose the properties you want for the state-space model. I would encourage you to read the Example10 notebook on the FreeGSNKE repo to see the different settings below. It basically determines how many passive structure modes you wish to use in the model - below I've selected to use all the modes (which is a bit slow). " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from freegsnke import nonlinear_solve\n", "\n", "state_space_model = nonlinear_solve.nl_solver(\n", " eq=eq, \n", " profiles=profiles,\n", " GSStaticSolver=GSStaticSolver,\n", " plasma_resistivity=2e-7, # set an appropriate plasma resistivity (goes into the Ip equation)\n", " # min_dIy_dI=1e-5, # set between [0,1]\n", " # threshold_dIy_dI=1, # set between [0,1]\n", " # max_mode_frequency=1e10, # set between [0,inf]\n", " fix_n_vessel_modes=138, # retain all passive modes in the MAST-U-like machine\n", " # custom_coil_resist=coil_R, # use calibrated resistance matrix for active coils\n", " # custom_self_ind=coil_M, # use calibrated inductance matrix for active coils\n", " n_linearization_workers=6, # set >1 to parallelise Jacobian construction (i.e. number of cores to use)\n", " )" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print(f\"Deformable GR = {state_space_model.linearised_sol.growth_rates}\")\n", "print(f\"Deformable GR (without Ip) = {state_space_model.linearised_sol.growth_rates_const_Ip}\")\n", "print(f\"Inductive stability margin = {state_space_model.linearised_sol.stability_margin}\")\n", "print(f\"Leuer parameter = {state_space_model.Leuer_metals_stab_over_active_destab}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now let's extract the matrices used in the FreeGSNKE state space model. Below is a (lengthy) recap of the equations that we use:\n", "\n", "\n", "---\n", "\n", "Assume we have solved the Grad-Shafranov equation and have obtained an equilibrium. From this we can extract the plasma current density discretised over the spatial domain (whose points are only non-zero inside the LCFS, so we mask points outside the limiter to save memory):\n", "\n", "\\begin{equation*}\n", " \\vec{I}_y := \\vec{I}_y^{GS} (\\vec{I}_m, I_p, \\vec{\\theta}).\n", " \\tag{1}\n", "\\end{equation*}\n", "\n", "This is a function of the: \n", "- currents in the metal conductors $\\vec{I}_m$ (i.e. the active coils and passive - or passive modes if we've done the decomposition).\n", "- plasma current $I_p$.\n", "- parameters $\\vec{\\theta}$ used to define the plasma current density profiles (i.e. $p'(\\psi_n)$ and $FF'(\\psi_n)$). \n", "\n", "In order to estimate the instability timescales/growth rates of this plasma, we will need to define (and then linearise) the coupled **circuit** and **plasma** equations. \n", "\n", "**Circuit equations**: \n", "\n", "These equations govern the flow of current in the active coils and passive structures (while being coupled to the currents in the plasma):\n", "\n", "\\begin{equation}\n", " M_{m,m} \\vec{\\dot{I}}_m + M_{m,y} \\vec{\\dot{I}}_y + R_{m,m} \\vec{I}_m = \\vec{V}_m\n", " \\tag{2}\n", "\\end{equation}\n", "\n", "- $M_{m,m}$ is the symmetric matrix of mutual inductances between all of the metals. \n", "- $M_{m,y}$ is the matrix of mutual inductances between all of the metals and the discretised plasma current density. \n", "- $R_{m,m}$ is the diagonal matrix of resistances in the metals.\n", "- $V_m$ is the vector of voltages applied to the metals (note: these are non-zero for the active coils only).\n", "- Note that the dots indicate time derivatives. \n", "\n", "**Plasma equation**: \n", "\n", "This equation governs the flow of current in the plasma (while being coupled to the currents in the metals):\n", "\n", "\\begin{equation}\n", " \\frac{\\vec{I}_y^T}{I_p} \\left( M_{y,y} \\vec{\\dot{I}}_y + M_{y,m} \\vec{\\dot{I}}_m + R_{y,y} \\vec{I}_y \\right) = \\vec{0}\n", " \\tag{3}\n", "\\end{equation}\n", "\n", "- $M_{y,m} = M_{m,y}^{T}$ (from above).\n", "- $R_{y,y} = 2 \\pi \\sigma_p \\hat{R} /dA$ is the diagonal matrix of resistances for each plasma element.\n", " - $\\sigma_p$ is the assigned plasma resistivity. \n", " - $\\hat{R}$ is a matrix storing the radial position of each plasma element.\n", " - $dA$ is the area of each plasma element.\n", "\n", "\n", "**Linearisation of discretised plasma current density**:\n", "\n", "To simplify these equations, we linearise equation (1) around the current GS equilibrium such that:\n", "\n", "\\begin{equation}\n", " \\dot{\\vec{I}_y} \\approx \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} \\dot{\\vec{I}}_m + \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS} \\dot{I}_p + \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}}.\n", " \\tag{4}\n", "\\end{equation}\n", "\n", "The partial derivative terms (i.e. the Jacobians) are calculated in FreeGSNKE using finite differences. The use of these Jacobians is what makes this a **deformable** growth rate model. \n", "\n", "We can plug this equation into both (2) and (3) and re-arrange:\n", "\n", "\\begin{align}\n", " \\left( M_{m,m} + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} \\right) \\vec{\\dot{I}}_m + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS} \\dot{I_p} + R_{m,m} \\vec{I}_m + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}} = \\vec{V}_m \\tag{5} \\\\\n", " \\frac{\\vec{I}_y}{I_p} \\left[ \\left( M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} + M_{y,m} \\right) \\dot{\\vec{I}}_m + M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS} \\dot{I_p} + R_{y,y} \\vec{I}_y + M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}} \\right] = \\vec{0}. \\tag{6}\n", "\\end{align}\n", "\n", "Next, we need to write this in state matrix form. For this, we need to remove the dependence of equation (6) on the denominator $I_p$ by introducing two terms:\n", "- $\\hat{\\vec{I}}_y = \\vec{I}_y / I_p$ is the normalised discretised plasma current density over the spatial domain.\n", "- $R_p = \\hat{\\vec{I}}_y^{T} R_{y,y} \\hat{\\vec{I}}_y$ is the \"lumped\" plasma resistance.\n", "\n", "This results in the state matrix system $\\vec{M} \\dot{\\vec{x}} + \\vec{R} \\vec{x} = \\vec{c}$:\n", "\n", "\\begin{equation*}\n", "\\begin{aligned}\n", "\\left[\\begin{array}{cc}\n", " M_{m,m} + M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} & M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS} \\\\ \n", "\\frac{\\hat{\\vec{I}}_y}{R_p} \\left( M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{I}_m} \\Bigg|_{GS} + M_{y,m} \\right) & \\frac{\\hat{\\vec{I}}_y}{R_p} M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial I_p}\\Bigg|_{GS}\n", "\\end{array}\\right]\n", "\\left[\\begin{array}{c}\n", "\\dot{\\vec{I}}_m \\\\ \n", "\\dot{I}_p\n", "\\end{array}\\right] +\n", "\n", "\\left[\\begin{array}{cc}\n", "R_{m,m} & 0 \\\\ \n", "0 & 1\n", "\\end{array}\\right]\n", "\\left[\\begin{array}{c}\n", "\\vec{I}_m \\\\ \n", "I_p\n", "\\end{array}\\right]\n", "\n", "&=\n", "\\left[\\begin{array}{c}\n", "\\vec{V}_m \\\\ \n", "0 \n", "\\end{array}\\right]\n", "-\n", "\\left[\\begin{array}{c}\n", "M_{m,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}} \\\\ \n", "\\frac{\\hat{\\vec{I}}_y}{R_p} M_{y,y} \\frac{\\partial \\vec{I}_y}{\\partial \\vec{\\theta}}\\Bigg|_{GS} \\dot{\\vec{\\theta}}\n", "\\end{array}\\right].\n", "\n", "\\end{aligned}\n", "\\end{equation*}\n", "\n", "Note that in the next section we use $\\vec{I}_d$ to represent the vector of active coil and passive *mode* currents such that $\\vec{I}_m = P \\vec{I}_d$ (where $P$ is the eigenmode decomposition matrix - see original freegsnke paper for details or Example10 notebook). \n", "\n", "Also note that you can ignore the effect of the final term on the right hand side of the equation (involving Jacobians wrt the plasma parameters) if so desired - these effects are thought to be small. \n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**State space model**:\n", "\n", "This can be written as:\n", "\n", "$$ L_* \\frac{d \\vec{x}}{dt} + R \\vec{x} = \\vec{u} $$\n", "\n", "where:\n", "- $\\vec{x}$ is the state space vector. \n", "- $L_*$ is the modified inductance matrix.\n", "- $R$ is the resistance matrix. \n", "- $u$ is the control (voltage) vector. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's look at the state space vector, which is made up of the:\n", "- 12 active coil currents [A]. \n", "- 138 passive structure currents [A] (assuming all modes being used). \n", "- plasma current [kA] (single number).\n", "\n", "The exact ordering is available from the machine object below." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# state space vector\n", "x = state_space_model.currents_vec\n", "print(x.shape)\n", "\n", "# ordering of the actives, passives, and Ip in the vector\n", "ordering = tokamak.coil_names\n", "ordering" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# mutual inductances between all metals\n", "Mmm = state_space_model.evol_metal_curr.coil_self_ind\n", "print(Mmm.shape)\n", "\n", "# resistances for all metals\n", "Rmm = state_space_model.evol_metal_curr.coil_resist\n", "print(Rmm.shape)\n", "\n", "# mutual inductances between all metals and plasma grid elements\n", "Mmy = state_space_model.evol_metal_curr.Mey_matrix\n", "print(Mmy.shape)\n", "\n", "# normalised plasma current density (Iy / Ip), sums to 1, defined over the same plasma-grid elements as Iy\n", "hatIy = state_space_model.hatIy\n", "print(hatIy.shape)\n", "\n", "# Myy @ hatIy, precomputed by FreeGSNKE without ever forming the full Myy matrix\n", "# (Myy itself is only ever built/stored on a small adaptive sub-domain around the plasma - see state_space_model.handleMyy)\n", "Myy_hatIy = state_space_model.Myy_hatIy0\n", "print(Myy_hatIy.shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# mode transformation matrix (if number of passive modes reduced)\n", "P = state_space_model.linearised_sol.P\n", "P.shape" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Jacobian of plasma current density wrt mode currents (i.e. the linearisation!)\n", "dIydId = state_space_model.linearised_sol.dIydI[:,0:-1]\n", "print(dIydId.shape)\n", "\n", "# Jacobian of plasma current density wrt plasma current (i.e. the linearisation!)\n", "dIydIp = state_space_model.linearised_sol.dIydI[:,-1] / state_space_model.plasma_norm_factor\n", "print(dIydIp.shape)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# diagonal of R_yy, the plasma resistance matrix (2*pi*sigma_p*Rhat/dA), one entry per plasma grid point\n", "Ryy = state_space_model.plasma_resistance_1d\n", "print(Ryy.shape)\n", "\n", "# lumped plasma resistance: Rp = hatIy^T @ Ryy @ hatIy\n", "Rp = np.dot(Ryy, hatIy**2)\n", "print(Rp)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Assemble $L_*$ and $R$\n", "\n", "Substituting $\\vec{I}_m = P \\vec{I}_d$ into the block matrix above (so that the state vector becomes $\\vec{x} = (\\vec{I}_d, I_p)^T$, matching `state_space_model.currents_vec`) gives:\n", "\n", "$$\n", "L_* = \\left[\\begin{array}{cc}\n", "M_{m,m}P + M_{m,y}A & M_{m,y}B \\\\\n", "\\dfrac{1}{R_p}\\left(\\hat{\\vec{I}}_y^{T} M_{y,y} A + \\hat{\\vec{I}}_y^{T} M_{y,m}P\\right) & \\dfrac{1}{R_p}\\hat{\\vec{I}}_y^{T} M_{y,y} B\n", "\\end{array}\\right],\n", "\\qquad\n", "R = \\left[\\begin{array}{cc} R_{m,m}P & 0 \\\\ 0 & 1 \\end{array}\\right]\n", "$$\n", "\n", "where $A = \\partial \\vec{I}_y/\\partial \\vec{I}_d$ (`dIydId`) and $B = \\partial \\vec{I}_y/\\partial I_p$ (`dIydIp`).\n", "\n", "Two things to watch for when using the quantities already extracted above:\n", "\n", "- $\\hat{\\vec I}_y^{T} M_{y,y} = (M_{y,y}\\hat{\\vec I}_y)^{T}$ since $M_{y,y}$ is symmetric, so both bottom-row terms can be built from the vector `Myy_hatIy0` alone (`Myy_hatIy0 @ A` and `Myy_hatIy0 @ B`) - we never need the full $M_{y,y}$ matrix." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "n = P.shape[0] # number of independent active-coil and passive-structure modes\n", "\n", "# L* blocks\n", "L_dd = Mmm @ P + Mmy @ dIydId # (n, n)\n", "L_dp = Mmy @ dIydIp # (n,)\n", "L_pd = (Myy_hatIy @ dIydId + (Mmy @ hatIy) @ P) / Rp # (n,)\n", "L_pp = (Myy_hatIy @ dIydIp) / Rp # scalar\n", "\n", "L_star = np.zeros((n + 1, n + 1))\n", "L_star[:n, :n] = L_dd\n", "L_star[:n, -1] = L_dp\n", "L_star[-1, :n] = L_pd\n", "L_star[-1, -1] = L_pp\n", "\n", "# R blocks\n", "R = np.zeros((n + 1, n + 1))\n", "R[:n, :n] = np.diag(Rmm) @ P\n", "R[-1, -1] = 1.0\n", "\n", "L_star.shape, R.shape" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Sanity check against FreeGSNKE's own growth rates\n", "\n", "$L_* \\dot{\\vec x} + R \\vec x = \\vec u$ implies $\\dot{\\vec x} = -L_*^{-1}R\\,\\vec x + L_*^{-1}\\vec u$, so the eigenvalues of $-L_*^{-1}R$ are the system's characteristic (growth) rates. \n", "\n", "This should match `state_space_model.linearised_sol.growth_rates` (up to which of the eigenvalues are unstable/positive), independently of the internal basis/scaling FreeGSNKE itself uses to build `linearised_sol.Mmatrix`.\n", "\n", "Note that FreeGSNKE internally calculates the GRs using $-R^{-1} L_*$ as it's more numerically stable so we'll do that here. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# P inverse\n", "Pm1 = state_space_model.evol_metal_curr.Pm1\n", "\n", "# find R^{-1} @ L* \n", "Rinv_L = np.zeros_like(L_star)\n", "Rinv_L[:n, :] = (Pm1 @ np.diag(1.0 / Rmm)) @ L_star[:n, :] # R_dd^-1 @ [L_dd | L_dp]\n", "Rinv_L[-1, :] = L_star[-1, :] # R_pp = 1, so this row is unchanged\n", "\n", "# eigenvalues of R^-1 L* are reciprocals of eigenvalues of L*^-1 R, so growth_rate = -1/eigenvalue\n", "timescales = -np.linalg.eigvals(Rinv_L)\n", "growth_rates = np.sort(1.0 / timescales[timescales.real > 0].real)\n", "\n", "print(\"Growth rates from L_*, R assembled here:\", growth_rates)\n", "print(\"Growth rates from state_space_model: \", state_space_model.linearised_sol.growth_rates.real)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "freegsnke", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.10.20" } }, "nbformat": 4, "nbformat_minor": 4 }